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ABSTRACT 

The  Burgers  equation  with  a  small  viscosity  term,  initial  and  periodic  boundary  condi¬ 
tions  is  resolved  numerically  using  a  spatial  approximation  constructed  from  an  orthonormal 
basis  of  wavelets. 

The  algorithm  is  directly  derived  from  the  notions  of  multiresolution  analysis  and  tree 
algorithms.  Before  the  numerical  algorithm  is  described  these  notions  are  first  recalled. 
The  method  uses  extensively  the  localization  properties  of  the  wavelets  in  the  physical  and 
Fourier  spaces.  Moreover,  we  take  advantage  of  the  fact  that  the  involved  linear  operators 
have  constant  coefficients.  Finally,  the  algorithm  can  be  considered  as  a  time  marching 
version  of  the  tree  algorithm. 

The  most  important  point  is  that  an  adaptative  version  of  the  algorithm  exists:  it  allows 
one  to  reduce  in  a  significant  way  the  number  of  degrees  of  freedom  required  for  a  good 
computation  of  the  solution. 

Numerical  results  and  description  of  the  different  elements  of  the  algorithm  are  provided 
in  combination  with  different  mathematical  comments  on  the  method  and  some  comparison 
with  more  classical  numerical  algorithms. 
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1.  INTRODUCTION 


The  nonlinear  parabolic  equation 


(1) 


du  udu  vd2u 
dt  dx  dx2 


known  as  the  Burgers  equation  is  one  of  the  simplest  combining  both  nonlinear  propagation 
and  diffusive  effects.  It  represents  a  first  step  in  the  hierarchy  of  approximations  of  the 
Navier-Stokes  equation.  Solutions  of  this  equation  exhibit  a  delicate  balance  between  the 
nonlinear  advection  and  the  diffusion  terms.  This  situation  leads  to  solutions  with  rapid 
and  localized  variations.  Moreover,  exact  solutions  are  known  (thanks  to  the  Hopf-Cole 
transformation).  Thus,  this  equation  appears  to  be  a  very  convenient  and  useful  test  problem 
for  new  numerical  schemes. 

This  paper  is  devoted  to  the  numerical  resolution  of  the  Burgers  equation  (1)  with  pe¬ 
riodic  boundary  conditions  on  [0,1]  and  known  initial  condition:  u(0,t)  =  u(l,  t)\ u(x,  0)  = 
sin(27rx).  We  use  a  new  algorithm  based  on  a  spatial  approximation  provided  by  an  orthonor¬ 
mal  wavelet  basis  exploiting  as  much  as  possible  the  numerical  localization  and  regularity 
properties  of  the  wavelets.  One  of  the  greatest  advantages  of  this  algorithm  is  that  the 
whole  method  can  automatically  adapt  itself  to  the  generation  of  large  gradient  regions  in 
the  solution. 

After  a  short  review  of  the  notions  of  multiresolution  analysis  and  tree  algorithms  that 
are  required  for  the  construction  of  orthonormal  wavelet  bases  and  that  are  the  keystones 
of  our  numerical  algorithm,  some  properties  of  r  —  regular  wavelet  bases  and  particularly 
of  spline  wavelets  (Section  2)  are  recalled.  The  choices  and  definitions  for  the  numerical 
algorithm  are  discussed  first  in  the  regular  case  (Section  3.1)  and  secondly  in  the  adapted 
case  (Section  3.2).  Numerical  results  are  presented  in  Section  4  with  a  complete  description 
of  the  different  elements  of  the  method.  Comparison  with  a  classical  Fourier  pseudospectral 
method  is  provided.  Section  5  is  devoted  to  concluding  remarks. 


2.  CONSTRUCTION  AND  BASIC  PROPERTIES  OF  ORTHONORMAL 
WAVELET  BASIS 


The  natural  framework  in  which  to  construct  bases  of  wavelets  is  given  by  the  multireso¬ 
lution  analysis.  This  new  concept,  elaborated  by  S.  Mallat  and  Y.  Meyer,  is  also  well  adapted 
to  the  approximation  of  functions  for  numerical  purposes. 

The  existence  of  a  “Fast  Wavelet  Transform,”  based  on  the  class  of  tree  algorithms  studied 
by  S.  Mallat  and  I.  Daubeschies,  is  also  a  nodal  point. 
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We  briefly  describe  the  main  definitions  and  results  we  will  need.  The  reader  is  referred 
to  the  book  of  Y.  Meyer  [Meyer  1987]  for  a  complete  account  on  the  subject.  Local  analysis 
properties  will  demonstrate  the  quality  of  the  wavelet  approximation  and  help  to  explain 
the  flagging  process  used  in  the  adaptative  version  of  the  algorithm.  The  periodic  spline 
wavelets  will  be  presented  at  the  end  of  this  section  and  the  different  properties  motivating 
their  choice  for  the  numerical  resolution  of  partial  differential  equations  will  be  recalled. 

2.1.  Multiresolution  Analysis 


As  introduced  by  S.  Mallat  and  Y.  Meyer,  a  multiresolution  of  L2(Mn)  is  a  specific 
approximation  scheme  for  finite  energy  functions.  In  this  paper  we  restrict  ourselves  to 
n  =  1.  A  multiresolution  of  L2(1R)  is  an  increasing  sequence  of  closed  linear  subspaces 
Vj,j  E  2Z  such  that: 

+oo  +oo 

(2.1)  =  ik  is  dense  in  L2(JR ) 

—  oo  —  oo 

(2.2)  f(x)  €  V,  <=►  /( 2x)  6  Vjtl 

(2.3)  fix)  EV.ts.  fix  -  k)  €  V0,  Vi  £  7Z, 

(2.4)  There  exists  a  function  g{x )  in  V0  such  that  the  collection  g(x  —  k),k  E  2Z  is  a  Riez 
basis  for  Vo 

(A  collection  of  vectors  ej,j  &  J  ,  in  a  Hilbert  space  H  is  a  Riez  basis  if  any  vector  x  E  H 
can  be  written  in  a  unique  way  as  a  sum  x  =  ej  where  (X)  lojl2)1^2  is  finite  and  defines 
an  equivalent  norm.) 

The  multiresolution  analysis  is  called  r-regular  if  the  function  g(x)  defined  by  (2.4)  has 
the  following  property: 

(2.5)  \dag(x)\  <  Cp{  1  +  |x|)_p  for  a  <  r,  all  x  E  R  and  all  integer  p. 

This  property  can  be  interpreted  as  a  quantification  of  the  localization  properties  of  the 
function  g  in  the  physical  space  (a  =  0)  and  in  the  Fourier  space  (a  >  1). 

Let  us  give  two  examples  of  multiresolution  analysis  that  generate  classical  spaces.  The 
first  one,  to  which  people  familiar  with  spectral  methods  automatically  think,  is  provided  by 
the  Paley-Wiener  analysis.  In  this  analysi  s,  Vj  is  made  of  functions  whose  spectrum  belongs 
to  the  interval  [— 2-7-1, 2J-1].  g  is  the  cardinal  sinus  g(x)  =  *in^x) .  It  is  well-known  that  g 
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is  infinitely  derivable  but  is  badly  localized  in  the  physical  space.  Consequently,  for  every 
value  of  r,  this  multi-resolution  analysis  is  not  r-regular. 

The  second  example  is  the  spline  multiresolution  which  will  provide  the  spline  wavelets 
that  will  be  used  in  the  following  numerical  implementation  of  the  algorithm  .  Here,  given  an 
integer  m  >  2,  Vj  consists  of  all  functions  with  m— 2  continuous  derivatives  whose  restrictions 
to  any  interval  [A:2— J ,  ( k  +  1)2-J]  coincide  with  a  polynomial  of  degree  less  or  equal  tom-1. 
g(x)  is  given  by  the  so-called  basic  spline  gm( x)  which  is  the  m  fold  convolution  of  the 
characteristic  function  of  [0, 1].  gm  is  supported  by  [0,m]  and,  satisfies  (2.5)  for  r  =  m  —  2. 

2.2.  The  Construction  of  Wavelets  from  the  Multiresolution  Analysis 

The  construction  of  the  wavelet  basis  stems  from  the  fact  that  during  the  process  of 
refinement  in  the  approximation  one  wants  to  only  store  the  improvement  from  the  approx¬ 
imation  j  to  the  approximation  j  +  1.  Mathematically,  one  introduces  at  each  step  j,  the 
subspace  Wj,  defined  as  the  orthogonal  complement  of  Vj  in  Vj+i.  The  Wj  space  family 
satisfies  the  scaling  (2.2)  and  translation  invariance  (2.3)  properties  imposed  on  Vj,  so  that 
attention  can  be  focused  on  W0.  Then  one  has  the  fundamental  theorem  proved  by  S.  Mallat 
and  Y.  Meyer. 

Theorem:  There  exists  a  function  xp  of  Wo  such  that  {xp(x  —  k),k  E  ZZ\  is  an  orthonormal 
basis  ofW0.  ip  has  the  same  regularity  properties  (2.5)  as  g,  and  xp  is  an  oscillating  function: 

(2.6)  J  xkxp(x)dx  =  0  for  0  <  k  <  r. 

Then,  the  family  of  functions  \xpjk{x)  =  2^2xp(2-’x  —  k),k  E  2Z}  is  an  orthonormal  family 
of  Wj  and,  as  L2(1R )  =  (BjezWj,  {ipjk(x),j  E  Z,k  E  Z}  is  a  hilbertian  basis  of  L2(M). 
This  family  is  an  orthonormal  wavelet  family.  Combined  with  (2.5)  for  a  >  0,  equation  (2.6) 
characterizes  the  localisation  of  the  wavelet  generating  function  in  the  Fourier  space. 


The  proof  of  the  theorem  introduces  an  orthonormal  basis  of  Vj  that  will  be  denoted 
by  {<pjk(x)  =  2*l2<p{2*x  —  k),  k  E  2S}.  These  functions  are  called  by  S.  Mallat  the  scaling 
functions  due  to  the  fact  that  /  <p(x)dx  =  1.  They  satisfy  the  same  regularity  properties  (2.5) 


as  g. 


Given  an  integer  j0  and  writing  L2(1R )  =  V]a  ©  Wy  one  obtains  another  hilbertian 
basis  of  L2{JR)  :  {<pj0,k1'Pjk,j  >  jo,k  E  Z). 

Let  us  go  back  now  to  the  two  examples: 

The  Paley-Wiener  Analysis  provide^  the  wavelet  generating  function  i p{x)  given  by: 


xp{x)  = 


2  cos  7 rx  —  sin  27rx 
7T  2x  —  1 


and  ip(w)  =  e  "rh'xi/2<|u-|<i 
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where  ip  stands  for  the  Fourier  transform  of  ip  :  ip( w )  =  /  ip{x)e~2l™xdx. 

A 

ip  has  a  compact  support  but  the  localization  of  ip  is  poor  as  it  decreases  as  ^  at  infinity 
and  ip(w)  =  0  if  |w|  <  |,  which  means,  in  a  weak  sense,  that  all  the  moments  of  ip  are 
vanishing. 

The  spline  multiresolution  analysis  provides  the  “exponentially  localized  wavelets”  de¬ 
rived  by  P.  G.  Lemarie  and  G.  Battle  [Battle  1987,  Lemarie  1989].  The  generating  wavelet 
ip(x)  satisfies  the  following  properties: 

(2.7)  ’P(x)  and  all  its  derivatives  to  the  order  m  —  2  are  continuous 


(2.8) 


there  exist  e  >  0  such  that  |-0(a:) |  <  Cce-*'^,  |^(x)|  <  Cie  *'*l,  •  •  ■ , 

OX 


dm~ 2tP{x) 
dxm~2 


<  Cm_9e-* 


-1  » 


sm 


/+°o 

x*ip(x)dx  =  0  for  k  =  0,  •  •  • ,  m  —  1 

-OO 

A 

ip  is  given  in  the  Fourier  space  by 

jr  s  e~iynt  [  Pm^i(coB8^) 

(^)m  (sin2  (^r))  Pm- 1  (sin2(7ritf)) 

where  Pm  is  the  m  order  polynomial  given  by 

Pm- i (sin2  to)  _  ^  1 

(sinu;)2m  ( w  +  A:7r)2m  ’ 

The  scaling  function  <p  is  also  “exponentially  localized”  and  satisfies: 


2m 


/  TTW\ 

\T) 


4>(w)  = 


l 


sin 


l(7 ru>) 


{kw)™  [Pm_1(sin2(7ru;))]i 
For  m  =  6  the  generating  spline  wavelet  ip  and  scaling  function  <p  are  plotted  in  Figure  1. 


2.3.  Tree  Algorithms 

The  tree  or  pyramidal  algorithms  are  the  basic  tools  for  the  fast  computation  of  wavelet 
coefficients.  Also,  the  numerical  algorithm  presented  here  is  efficient  thanks  to  its  pyramidal 
structure. 
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If  Jm  is  a  given  integer,  fixing  the  level  of  approximation  at  2_JM,  the  tree  algorithm 
allows  fast  computation  of  the  wavelet  coefficients  from  jm  —  1  to  an  order  jo,  Jo  <  jM, 
starting  from  the  scaling  coefficients  in  VjM. 

Writing 

1 

(2-10)  Viu  =  VJ,  ®  £  W, 

i-jo 

and  starting  from  the  j^  approximation  of  a  function  U,  we  have 

UvjM{U){x)  =  Y,  2,m/2cJm (j){VMx  -  k). 
kez 

(Generally,  Ft*  stands  for  the  orthogonal  projector  from  L2(M )  on  the  subspace  X.)  At 
each  level  j,  S.  Mallat  [Mallat  1988]  shows  that,  due  to  the  orthogonality  properties  and  the 
scaling,  IIiyi_l(C/')  and  II ^.^(17),  written  as: 

nwi_1(f7)(a:)  =  53  2 ^  *j-xjfll>{V~Xx  -  i) 

l€Z 

and 

tV,(C0W  =  £  2i^ci_,^(2>-1i  -t) 

tez 

can  be  computed  from  11^(1/)  using  the  formulas: 


dj-1,1  —  cj<k  dfit  k) 

k 


—  53  cj,k  h(2i  k ) 

k 

where  g  and  h  are  two  discrete  functions  depending  only  on  the  multiresolution  analysis  and 
independent  of  j.  These  functions  are  precisely  defined  as: 

Vn  e  N  :  g(n)  =  -J=  <  V’-i.o,  <t>o,n  >  and  h[n )  =  -^=  <  4>- 1,0>  <Po,n  > 


9{~n)  =  9(n)  ,  h(-n)  =  h(n) 


where  <  ., .  >  stands  for  the  L2{IR)  inner  product  (<  f,9  >=  J  f(x)g(x)dx). 

Practically,  it  means  that  with  the  storage  of  two  discrete  functions  h  and  g  which  are  well 
localized  thanks  to  the  localization  of  ip  and  <p  ,  one  computes  with  an  order  of  iVLogjV(2,M  = 
N)  operations  the  (2;M— 2;o)  wavelet  coefficients  and  the  2:o  scaling  coefficients  corresponding 
to  the  formula  (2.10).  The  discrete  functions  h(n)  and  g(n)  are  plotted  on  Figure  2  in  the 
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physical  and  Fourier  spaces  in  the  case  of  the  spline  multiresolution  with  m  =  6.  They 
classically  appear  as  c  low  pass  filter  and  a  high  pass  filter  respectively. 

2.4.  Approximation  properties  of  Wavelets 

Approximation  properties  of  orthonormal  wavelets  can  be  interpreted  from  different 
points  of  view: 

As  the  space  generated  by  the  wavelets  up  to  an  order  j  —  1  is  the  space  Vj  of  the 
multiresolution,  the  approximation  results  are  driven  by  the  properties  of  the  scaling  function 
2,/a^(2,x)  that,  with  its  translates  <f>jk ,  generate  Vj.  It  is  then  a  more  classical  problem  of 
approximation  by  translates.  The  general  results  relating  the  quality  of  the  approximation 
to  the  degree  r  of  r-regularity  of  the  multiresolution  can  be  found  in  Y.  Meyer’s  book  [Meyer 
1990]. 

Moreover,  if  the  space  Vj  is  split  in  the  different  contributions  of  Wi,  l  <  j  —  1,  a  char¬ 
acterization  of  the  Sobolev  space  Ha(M )  can  be  obtained  from  the  decay  of  the  wavelet 
coefficients  as  : 

Theorem:  If /belongs  to  H~T(M),  ifVj,j  £  2Z  is  a  r-multiresolution  of  L2(IR)  and  if  — r  < 
s  <  r  then  fbelongsto  H'(M)  if  and  only  i/IIy0(/)  €  L2(1R )  and  |  |IItv,(/)|  b  =  j  €  IN, 

where  £j  €  l2  (IN). 

Another  important  result  for  approximation  by  wavelets  deals  with  the  characterization 
of  singularities  from  the  wavelet  decomposition.  The  singularity  of  a  function  at  a  point  Xq 
can  be  described  with  the  Lipschitz  exponent  a  defined  as: 

f(x)  is  a  Lipschitz  at  xo  ,(0  <  a  <  1)  if  and  only  if  for  all  x  in  a  neighborhood  of  xo  one 
has: 

l/(x)  —  f(xo)\  =  O(\\x-x0\n 
The  following  theorem  holds: 

Theorem:  A  function  f  belonging  to  L2(M )  is  a  Lipschitz  with  0  <  a  <  1  in  all  points  of 
an  open  interval  if  and  only  if  for  all  x  in  the  interval  one  has: 

I  </,*#>  I  =  0(2-' 


As  pointed  out  by  Stephane  Mallat  [Mallat  1990],  an  extension  of  this  result  can  be 
obtained  when  the  initially  a  singular  function  f  is  smoothed  by  a  smoothing  kernel.  For 
instance  with  a  gauusian  of  variance  a  one  has  the  following  result: 
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The  wavelet  coefficient  of  the  smoothed  function  at  the  scale  2~3  is  of  the  same  order  of 
magnitude  as  the  wavelet  coefficient  of  the  original  function  f  at  the  scale  s  =  \/2~23  +  a2. 
The  previous  theorem  then  holds  by  changing  2~ 3  in  \J2~ 23  +  a2. 

On  one  hand,  when  2~3  o,  i.e.,  when  the  scales  under  study  are  larger  than  the 
smoothing  scale,  the  previous  theorem  shows  that  the  characterization  of  the  degree  a  from 
the  wavelet  coefficients  of  the  smoothed  function  is  still  possible.  On  the  other  hand,  when 
2~3  <  a,  the  scaling  effect  disappears  and  no  characterization  of  the  degree  of  singularity  is 
possible  from  the  wavelet  coefficients  of  the  smoothed  function. 

2.5.  Periodic  Spline  Wavelets 

As  pointed  out  by  Y.  Meyer  [Meyer  1990],  the  complete  tool  box  built  in  L2(M)  can  be 
used  in  the  periodic  case  for  L2([ 0, 1])  by  introducing  a  standard  periodization  technique.  It 
has  been  implemented  by  V.  Perrier  and  C.  Basdevant  [Perrier  and  Basdevant  1989].  This 
technique  consists  at  each  scale  in  folding,  around  the  integer  values,  the  wavelets  ipjk  and 
the  scaling  functions  <f>jk  centered  in  [0, 1]  resulting  in 

i>P:k(x)  =  52  -  l ) 

leZ 

and 

4>Pjk(x)  =  52  Mx  ~  0- 
KZ 

This  folding  provides  a  multiresolution  analysis  of  L2([0,1]).  Due  to  the  finite  size  of  [0,1], 
j  takes  only  positive  values  and  k  takes  a  finite  number  of  values  at  each  scale:  0  <  k  < 
23  —  1.  Moreover,  the  periodized  function  <ftpo,o  that  generates  Vpo  is  constant  and  equal  to 
1.  Following  V.  Perrier  and  C.  Basdevant,  the  periodic  m  order  spline  scaling  functions  and 
wavelets  are 

1  V’Pj.oH  =  3^  (fr) 

where  <frpjk  and  i’Pjk  generate  respectively  Vp:  and  Wpj. 

For  simplicity  in  the  notation,  from  now  on  the  subscript  p  will  be  removed.  When 
necessary  the  subscript  up  will  be  added  for  the  non  periodized  functions. 

As  far  as  numerical  application  is  concerned,  the  spline  multiresolution  offers  the  following 
advantages. 

First,  it  is  known  that  if  the  degree  m  of  the  spline  is  even,  there  exists  a  Lagrangian 
interpolant  function  gjM  of  VjM  such  that  gjM{k2~3M)  =  6k,o-  (6  stands  for  the  Kronecker 
symbol:  6k,o  =  1  if  and  only  if  k  =  0.)  The  functions  gjM(x  —  k2~3M )  0  <  k  <  23M  are  then 
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an  unconditional  basis  of  VjM.  Thanks  to  this  basis,  it  is  possible,  for  any  function  I/(x),  to 
consider  its  interpolation  on  VjM  given  by 

n.viM(£0(*)  =  £  U{k2-i«)g^{x-k2-’“). 

This  collocation  projection  is  very  useful  when  one  wants  to  apply  the  classical  but  efficient 

pseudospectral  treatment  to  nonlinear  terms  (see  Section  3.1.5). 

A  second  point  of  interest  stems  the  numerical  localization  of  the  generating  functions 

V’u p  and  (pup  in  physical  and  Fourier  spaces.  Numerically,  the  exponential  decay  of  these 

functions  can  be  described  by  the  required  support  to  which  ipup  or  <PuP  (or  xpup  and  <pup)  can 

be  restricted  without  altering  their  L2  norms  (equal  to  1)  up  to  a  fixed  precision  a. 

In  the  physical  space,  the  support  of  the  periodized  functions  <pjk  and  ipjk  is  directly 

connected  to  the  support  of  <pup  jk  and  ypuP  jk  as  soon  as  the  length  of  these  last  supports  is 

less  than  1  (otherwise  the  support  of  <pjk  and  ipjk  is  [0,1]).  Moreover  if  S^p  and  S^,up  are 

s  s 

the  supports  of  (f>up  and  ipup  the  supports  of  <pup  jk  and  i pup  jk  are  respectively  and 
For  m  =  6,  the  following  estimations  are  obtained: 


precision  a 

support  of  <PuP 

support  of  ipup 

IF3  1 

|x|  <  3.50 

|x |  <  3.25 

io-4 

jxj  <  4.75 

jxj  <  4.75 

and  an  estimate  of  the  increase  of  the  length  of  these  supports  with  m  can  be  found  in  [Perrier 
and  Basdevant  1989].  For  j  >  3,m  =  6  and  up  to  a  precision  a  =  10-4,  the  functions  (<pup)jk 
and  (V’up)jfc  have  a  numerically  compact  support  of  length  less  than  1. 

In  the  Fourier  space,  the  periodization  does  not  affect  the  support  of  the  functions  <p  and 
ip.  If  Si  are  the  bounds  of  the  support  of  p  (or  ip),  2* Si  ai^  the  bounds  of  the  support  of  <pjk 
(or  Tpjk)-  For  m  --  6  we  obtained 


precision  a 

support  of  (p 

support  of  -tp 

IF3  ™ 

□ 

tu|  <  0.625 

0.34  < 

w\  <  1.29 

10"4 

M  <  0.685 

0.27  < 

ltt/|  <  1.40 

The  length  of  these  supports  decreases  as  m  increases.  We  will  strongly  use  the  fact  that, 
for  m  =  6  and  a  =  10  ,  the  numerical  supports  of  ipjk  and  ipj+2,k  are  disconnected. 

A  third  advantage  of  the  spline  wavelets  is  the  very  easy  control  of  their  cancellation 
(the  first  m  —  2  moments  of  each  wavelet  vanish).  As  has  been  shown  in  Section  2.4,  on 
one  hand,  this  property  is  required  to  have  a  good  convergence  of  the  wavelet  series  and  a 
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good  characterization  of  the  local  properties  of  the  analyzed  function.  On  the  other  hand, 
as  m  increases  the  localization  in  the  physical  space  of  the  wavelet  decreases.  To  obtain  a 
“satisfactorily  localized  function”  with  “good  cancellation,”  one  has  to  define  a  compromise 
for  the  choice  of  m.  All  the  numerical  results  presented  in  this  paper  are  obtained  with 
m  =  6. 

To  be  complete,  one  must  recall  that  an  efficient  implementation  of  periodic  spline  wavelet 
decomposition  has  been  described  by  V.  Perrier  and  C.  Basdevant  [Perrier  et  al.  1989].  It 
allows  computation  with  order  of  NLog(N)  operations  of  the  wavelet  coefficients  of  a  periodic 
function  described  by  N  regularly  distributed  collocation  points.  This  algorithm  has  been 
extended  to  a  non  complete  (non  regular)  distribution  of  points. 

3.  THE  ALGORITHM 

Taking  apart  the  Lagrangian  methods,  numerical  algorithms  can  be  roughly  speaking 
separated  into  three  families:  finite  difference  schemes,  finite  element  methods  and  spectral 
schemes. 

The  class  of  finite  difference  schemes  is  known  to  be  very  efficient  in  intricate  configura¬ 
tions  because  it  can  be  easily  adapted  to  resolve  localized  difficulties  such  as  large  gradient 
regions  or  sharp  boundary  conditions.  This  flexibility,  however,  hides  a  difficulty  in  control¬ 
ling  precision.  Moreover,  the  computation  cost  increases  considerably  with  the  precision. 

The  finite  element  schemes  offer  a  better  controlled  precision  and  are  extensively  used 
for  intricate  boundary  condition  problems  in  multidimension.  However  they  suffer  from 
difficulties  in  getting  very  flexible  and  fast  algorithms. 

The  precision  of  spectral  methods  is  well-known  ([Gottlieb  and  Orszag  1977],  [Canuto 
et  al.  1987])  and  their  efficiency  for  problems  with  a  high  degree  of  regularity  has  been 
demonstrated.  However,  they  are  fundamentally  ill-suited  to  problems  that  develop  strong 
local  gradients  or  discontinuities  in  their  solutions. 

The  wavelet  approach  provides  a  multiscale  decomposition  based  on  orthonormal,  regular 
and  numerically  well  localized  functions.  It  can  then  a  priori  offer  an  interesting  compromise 
between  precision,  efficiency,  and  adaptability. 

3.1.  Description  of  the  Algorithm  on  a  Regular  Grid 


3.1.1.  The  Grid  Points 

Given  a  family  of  wavelets,  one  can  characterize  the  space  X  that  they  generate  by  the 
dyadic  grid  built  from  the  set  of  all  yjk  =  ^  +  ^rr  suc^  that  V'yfc  belongs  to  X  (see  Figure  3). 
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If  X  is  a  space  VjM,  then  a  regular  grid  corresponding  to  the  centers  of  the  scaling  functions 
4>juk  is  built  from  the  s*.  is  the  level  of  approximation. 

3.1.2.  The  Algorithm 


The  algorithm  is  first  described  generally  for  the  evolution  equation 


(3.1) 


'  %  +  LU  +  G(U)  =  0 

U(0,t)  =  U(l,t) 

< 

U  =  U0  for  t  =  0 
t  >  0,0  <  x  <  1 


where  L  is  a  linear  operator  and  G  is  a  function  of  U  and  its  derivatives.  An  approximate 
solution  Um,  belonging  to  a  trial  space  Xm,{Xm  =  VjM  in  the  case  of  a  regular  grid)  is 
sought  as  a  solution  of  the  equation: 

’  -af  +  LmUm  =  UXm(G(Um)) 

< 

.  uM(x,o)  --  njMUo. 

Here  LM  is  an  approximation  of  L,  of  the  form  PxMLTLxM ,  where  II^j,  is  the  orthogonal 
projection  onlo  XM,  and  PXm  is  another  projection  onto  Xm  parallel  to  the  orthogonal  of 
Ym,  another  space  oo  nx^CGCuM));  the  approximate  equation  becomes 


(3.2) 


+  LmUm  =  Gm 
UM(x,0)  =  U0m(x). 


Using  the  classical  name  conventions  for  the  Methods  of  Weighted  Residuals  (MWR)  (see 
[Canuto  et  al.  1987]),  we  call  trial  functions  the  wavelets  ipjk,0  <  j  <  <  k  <  2J  and 

<$o,o  and  introduce  a  family  of  functions  9t,t  €  T,  called  the  test  functions:  they  generate  our 
space  Ym .  The  weighted  residuals  minimization  statement  is  then  equivalent  to  the  following 
variational  formulation: 


(3.3) 


<  +  LmUm  —  Gm,v  >=  0  for  all  v  £Ym 

,  UM{x,0)  =  U0st(x) 


where  <  •,•  >  stands  for  the  L2  inner  product  on  [0,1].  In  the  case  of  the  Burgers  equation 
(1),  we  will  have  LU  =  G(U)  =  —  | ^ ,  Gm(U )  =  IIC(G({/)  and  the  test  functions 

6r  will  be  defined  in  the  next  paragraph. 
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3.1.3.  Choice  of  the  Test  Functions  and  of  the  Time  Discretization  Scheme 

The  general  idea  of  the  algorithm  is  to  comply  with  the  localization  of  the  wavelets  in 
the  Fourier  space  and  in  the  physical  space.  We  know  that,  in  comparison  to  the  Fourier 
functions  e,fcr,  the  localization  in  the  physical  space  implies  the  mixing  of  the  frequencies 
that  prevents  the  wavelets  from  being  eigenfunctions  of  differential  operators.  However,  they 
are  not  far  from  being  so.  Indeed,  the  localization  of  the  wavelets  in  the  Fourier  space  leaves 
the  differential  operators  nearly  diagonal. 

In  the  particular  case  of  the  spline  wavelets,  it  has  been  shown  previously  that,  to  a  given 
precision  a,  the  numerical  supports  of  V’jfc  and  ipjn  are  disconnected  as  soon  as  | j  —j'\  >  n(a) 
(for  a  =  10-3,n  =  1).  If  D  is  a  differential  operator  with  constant  coefficients,  the  same 
property  holds  also  for  DV'jfc  and  V’j'Jfc  with  a  new  precision  a1  depending  on  the  spectrum  of 
D.  If,  for  instance  D  =  (1  —  the  precision  is  unchanged. 

The  algorithm  we  are  going  to  describe  here  takes  large  advantage  of  this  property. 
Moreover,  thanks  to  the  fact  that  the  coefficients  of  the  Burgers  equation  are  constant,  some 
precalculations  can  be  performed  once  and  for  all.  Indeed,  the  main  part  of  the  work  to  be 
performed  at  each  time  step  of  the  MWR  can  be  done  only  one  time  and  reused  at  every 
time  step.  This  is  done  by  a  proper  choice  of  the  test  functions  6r. 

Let  us  be  more  precise  in  the  case  of  an  Euler  time  discretization  scheme,  implicit  for  the 
dissipation  term  and  explicit  for  the  convection  term.  This  choice  is  made  only  for  sake  of 
simplicity  in  the  following  presentation  and  more  sophisticated  schemes  are  currently  used 
for  the  numerical  implementations  (see  Section  4). 

Equation  3.3  is  discretized  as  follows: 

|  <  (1  -  ‘'te&W'.v  >=<  VJ,  -  >  for  all  v  S  Ym 

\  Uu(x,  0)  =  %„(*)• 

Choosing  r  =  (j,k)  and  0T  =  8jk  such  that 


(3.5) 


[  Mx)  =  {J  -  ^1?)  1  <Mx) 


[  80(x )  =  <f> 00(x) 

we  invert  once  and  for  all  the  above  direct  problem  at  each  scale. 
Following  this  choice,  the  following  set  of  equations  is  obtained: 

<  utt'A*  "=<  uz.Ojk  >  +  <  -f  0*  > 

<  ^M+\0OO  >  =  <  Um ,#0  > 
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that  can  be  written  as: 


f  <  utf'Ajk  >=<  >  +  <  (us,y,e'jt  > 

(3.6) 

i  <  ^Af+1>^00  >  =  <  Um,0o  > 
where  we  introduced  another  set  of  functions: 

***(*)  =  -f  (A  (7  “  vAt&))  1  V»jfc(®) 

(3.7) 

.  «o(*)  =  0. 

Using  a  Crank  Nicholson  discretization  scheme  and  spline  wavelets  of  order  m  =  6,  the 
function  6  associated  with  ^up  using  (3.5)  is  plotted  on  Figure  5.  Roughly  speaking,  the  Qjk 
functions  look  like  the  corresponding  wavelets  ipjk'.  In  the  case  of  the  spline  wavelets,  the 
Qjk  functions  are  exponentially  localized  in  both  Fourier  and  physical  spaces.  Moreover  they 
have  the  same  number  of  zero  moments  as  the  original  wavelets  and  their  numerical  support 
in  the  Fourier  space  is  included  in  those  of  the  corresponding  wavelets.  At  each  scale,  the 
Qjk  deduce  one  from  the  other  by  translation  but  the  wavelet  rescaling  property  from  one 
scale  to  the  other  is  no  longer  valid. 

3.1.4.  Computation  of  <  Qjk  > 

A  fast  computation  of  these  coefficients  is  possible  thanks  to  the  previous  remarks.  If, 
following  2.10,  one  writes: 

vj  <  j\i  ~  2  nVjM  —  nVj„  ©  53  nw, 

following  3.1.1  with  D  =  (1  —  i/A t-^),  a  numerically  good  approximation  of  <  U^,Qjk  >  is 
<  IIy,.+J(t/jk),  Qjk  >■  Then,  recalling  that: 

nvi+J  =  nVj+l  ©  uWj+l 

one  gets 

<  UhijQjk  >—<  Hv,+a(U^f), Qjk  >—<  nv,+i {Um)> Qjk  >  +  <  Uw,+1  Qjk  >  . 

Since  we  have  IIvj+l(U]^)  =  T,Cj+i,i<l>j+U  and  =  Erfj+i.W'j+U  we  obtain 

<  Um *  Qjk  53  Ci+M  “  2k )  +  dj+uPAZ  ~  2k) 

i 

where  a,  and  f3j  are  two  families  of  discrete  functions  parameterized  by  j  and  depending  on 
the  multiresolution  analysis,  on  the  equation  and  on  the  time  scheme.  They  are  defined  as 

Vn  €  IN  otj(n)  =<  <f>j,nfQj- i,o  >  and  / 3j(n )  =<  ipj,n,Qj- i,o  > 
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and 

a(— n)  =  aj(n )  and  /3j(—n )  =  /3j(n). 

As  in  the  tree  algorithm  (Section  2.3),  the  calculation  of  {<  Uj^,6jk  >,  for  0  <  j  < 
]M-\  and  0  <  Jfc  <  2>}  involves  two  filtering  procedures.  The  whole  calculation  requires 
order  JVLogJV  operations.  The  coefficients  of  the  niters  aj  and  /3j  are  stored  in  the  memory 
at  the  beginning  of  the  calculation  in  an  array  of  size  proportionnal  to  NLogN  . 

One  wants  to  emphasize  that  the  one-time  computation  involved  in  the  calculation  of  the 
functions  6jk  is  possible  because  the  operator  I  —  vAt-jj^  has  constant  coefficients. 

The  computation  of  <  Ufa,  Q'jk  >  can  be  done  the  same  way  using  two  other  families  of 
filters,  a'-,/?',  as  soon  as  (f/jj/)2  is  computed  (see  next  section). 

Recalling  that  the  filters  aj,  (3j,  a'j ,  /3'-  depend  on  At,  it  must  be  said  that  the  independence 
of  At  with  respect  to  j  is  not  required.  Practically  it  means  that  the  time  step  At  can  be 
adapted  to  the  spatial  scale,  a  possibility  that  can  be  useful  in  optimizing  the  time  stability 
and  precision  conditions. 

3.1.5.  Estimation  of  the  Nonlinear  Terms 


The  calculation  of  the  nonlinear  terms,  precisely  (Ufo)2,  is  done  using  the  classical  collo¬ 
cation  technique. 

3.1.6.  Implementation  of  the  Algorithm 

At  each  time  step,  the  algorithm  works  as  follows:  supposing  that  ( )  is  known  by  its 
values  at  grid  points,  its  scaling  coefficients  or  its  wavelets  coefficients,  U^1  is  computed 
using  an  order  of  NLogN  operations  as  shown  in  Table  1. 


UnMeV; 


Pseudospectral 

- ►  {UnM?  e  viM 


Tree  Algorithm 


Table  1. 
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3.2,  Description  of  the  Algorithm  on  a  Non  Regular  Grid 


Up  to  now,  Xm  was  equal  to  Vilt  which  means  that  every  wavelet  i/'jfc  with(j,  k )  such 
that  0  <  j  <  Jm  and  0  <  k  <  2J  —  1  was  considered  in  the  approximation  of  the  solution. 

However,  the  distribution  of  scales  significantly  building  a  given  function  in  the  sense 
provided  by  its  wavelet  decomposition  has  no  reason  to  be  regular  in  the  whole  interval 
[0,1].  Figure  7  shows  the  wavelet  decomposition  of  the  solution  of  (1)  as  large  gradient 
regions  develop  locally  around  x  =  .5.  According  to  the  values  of  the  wavelet  coefficients, 
small  scales  (j  >  5)  are  only  meaningful  in  a  region  very  close  to  .5  and  for  time  larger  than 
.5/ 7r.  More  precisely,  if  we  choose  a  given  precision  a  on  the  L2  norm  and  look  at  the  number 
of  wavelets  required  to  approximate  the  solution  at  different  times  it  appears  that  only  a 
small  number  of  coefficients  of  the  wavelet  decomposition  are  needed.  For  instance,  with  a 
L 2  norm  precision  of  10~6  the  number  of  wavelets  required  to  reconstruct  the  solution  at 
t  =  l/7r  is  74.  The  distribution  of  these  “active”  wavelets  follows  the  cone  shape  structure  of 
Figure  7.  Our  main  goal  is  to  compute  only  the  wavelet  coefficients  <  >  in  a  region 

as  close  as  possible  to  this  “active”  space. 

The  adaptative  version  of  our  algorithm  is  designed  in  order  to  fulfill  efficiently  this 
requirement.  It  involves  different  steps: 

Step  1:  Estimation  of  the  required  space  of  approximation  for  the  next  time  step:  In 
a  way  comparable  to  what  is  done  in  spectral  methods  (see  Gottlieb  -  Orzag  1977),  the 
accuracy  of  the  computation  is  examined,  locally  in  the  case  of  wavelets,  from  the  shape  of 
the  local  spectrum  provided  by  the  wavelet  decomposition  at  every  point. 

Step  2:  Adaptation  of  the  space  of  approximation  for  the  next  time  step:  According 
to  the  result  of  the  previous  estimate  and  the  knowledge  of  the  operators  of  the  Burgers 
equation,  the  space  of  approximation  Xm  is  optimized  (Xm0)  for  future  time  by  truncation 
of  non  “active”  small  scale  wavelets  or  local  addition  of  smaller  scale  wavelets  (that  are 
supposed  to  become  “active”  during  the  future  time  steps).  (The  smallest  scale  wavelets  in 
XMa  belong  to  WJMa_ i). 

Step  3:  Time  advancing  on  the  adapted  space  Xm*- 

We  are  going  now  to  describe  in  more  detail  the  adaptative  version  of  the  algorithm 
in  terms  of  flagging  criterion  (Step  1)  and  new  trial  and  test  functions  (Steps  2  and  3). 
Moreover  it  will  be  shown  that  the  general  structure  of  the  algorithm  remains  unchanged. 

3.2.1.  New  Trial  Space  Xm  and  Test  Space  Ym  but  Same  Algorithm 

Figure  4  shows  the  selected  values  of  ( j ,  k )  of  an  adapted  grid  generated  by  the  adapted 
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version  of  the  algorithm  for  the  Burgers  equation  (1)  at  t  —  1  /«.  Precisely, 

XMa  =  span  {ipjk,{j,  k)  G  JK} 
where 

JK  =  {( j ,  k )  such  that  0  <  j  <  j^,  0  <  k  <  2J 
or  such  that  j  =  jM,k\  <  k  <  k2 
or,  j  =  +  1,  k[  <  k  <  k'2) 

with  jm  =  6,  k\  —  22  k2  =  31  k[  =  54  k'2  =  73.  j^a  =  8. 

To  define  completely  the  algorithm,  one  has  also  to  define  the  Ym3  space.  This  is  done  in 
exactly  the  same  way  as  it  has  been  presented  in  the  regular  case:  the  generating  functions 
of  Yms  axe  derived  from  the  basis  functions  of  Xmo  using  (3.5).  Due  to  the  local  form  of 
the  algorithm,  the  whole  procedure  remains  unchanged  when  the  trial  and  test  spaces  are 
adapted.  The  complete  procedure  is  identically  reproduced  without  modification. 

This  highly  attractive  property  is  due  to  the  localization  properties  of  the  wavelets  and 
of  the  algorithm  in  physical  and  Fourier  spaces. 

The  adaptive  procedure  is  precisely  a  controlled  modification  of  the  space  of  approxima¬ 
tion  in  which  one  looks  for  the  solution.  As  this  procedure  consists  in  adding  or  subtract¬ 
ing  orthonormal  functions,  the  so-called  restriction  and  interpolation  operators  are  exactly 
known:  they  are  directly  provided  by  the  tree  algorithm.  Moreover,  no  reversibility  problems 
are  encountered  during  the  restriction  and  interpolation  procedures. 

For  completeness,  some  words  must  be  said  on  the  flagging  procedure  used  to  modify 
the  subscript  space  JK  (this  is  the  practical  way  to  adapt  Xm).  In  the  case  of  the  Burgers 
equation,  this  procedure  is  very  simple  and  cheap.  Basically,  starting  from  the  smallest  scale 
available  (jss)>the  attention  is  focused  on  the  smallest  scale  projection  of  the  solution 

n»i„(ua)  =  £  <  Vw  >  Vw- 

k  such  that  [jsSttyGJK 

As  shown  in  Section  2  and  Figure  8,  this  function  is  localized  in  the  strong  gradient 
regions.  The  energy  associated  with  the  wavelets  precisely  |  <  ipjSSlk,  Vh  >  I2, 

reveals  the  amount  of  energy  of  the  solution  in  the  scale  jss  around  the  point  yjjr*  + 
jjjz-  The  Burgers  equation  is  such  that,  due  to  the  presence  of  nonlinear  terms  in  the 
equation  that  generate  smaller  and  smaller  scale,  only  small  values  of  these  coefficients  are 
relevant  (otherwise  oscillations  occur).  Consequently,  if  |  <  t/>Jssfc,  Ufo  >  \ 2  is  greater  than 
a  predetermined  value,  the  space  Xm  of  approximation  is  adapted  by  adding  new  smaller 
scaled  wavelets  (ja  =  jss  +  1)>  around  the  point  5757^  +  5777 •  Adaptivity  means  also 
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restriction  of  the  space  XMa  and,  using  a  comparable  criterion,  small  scale  (jss)  wavelets  can 
be  removed  in  a  region  where  the  terms  |  <  Ufa,  ipjssi\2  are  not  significant.  (This  occurs  in  the 
second  regime  of  the  Burgers  solution  when  the  viscosity  smooths  the  gradients  previously 
generated  during  the  first  regime.) 

4.  NUMERICAL  RESULTS 

This  section  is  devoted  to  a  practical  and  numerical  description  of  the  algorithm.  The 
basic  elements  of  the  problem  and  the  algorithm  are  repeated  for  clarity. 

The  periodic  regularized  Burgers  equation 

'  §U  4.  U9U  _  i >d*U 
dt  '  dx  ~  dx3 

*>0,x<E[0,l],i/  =  ±^ 

i 

U(0,t)  =  U(l,t) 

U(x,0 )  =  sin(27rx) 

is  solved  numerically  in  the  space  of  splines  of  order  m  =  6.  The  trial  functions  are  the 
periodic  spline  wavelets  i/’iJt  and  <t> oo  •  The  test  functions  are  defined  following  the  general 
procedure  described  in  Section  3.1.5.  All  the  following  results  are  obtained  using  an  Adams- 
Bashforth  time  scheme  for  the  nonlinear  convection  term  and  a  Crank- Nicholson  time  scheme 
for  the  diffusion  term.  The  integration  time  step,  for  every  scale  is  At  —  10~3.  We  recall  here, 
that  the  choice  of  the  time  scheme  and  the  value  of  At  directly  influence  the  definition  of  the 
functions  ejk  and  6jk.  However,  knowing  the  time  discretization  scheme,  their  estimation  is 
8taightforward. 

The  filters  used  in  the  wavelet  decomposition  are  noted  h  and  g.  The  filters  defined  to 
compute  the  scalar  products  <  Uj^,  Ojk  >  and  <  (U^)2,  #'fc  >  are  noted  a j ,  /3j  and  /ir¬ 
respectively.  The  Figures  1,  2,  and  5  show  in  the  Fourier  and  in  the  physical  spaces  some 
selected  elements  of  the  different  mathematical  beings  defined  above. 

The  evolution  of  the  approximated  solution  U(x,t )  from  t  =  0  to  t  =  ^  is  plotted  on 
Figure  6  where  =  8  and  a  regular  grid  (Figure  3)  is  used.  Considering  the  theoretical 
solution  of  the  Burgers  equation  (see  [Basdevant  et  al.  1986]),  this  can  be  considered  as  a 
reference  estimation  of  the  solution.  Figure  7  shows  the  wavelet  coefficients  of  the  solution 
at  the  same  times.  Figure  8  shows  the  time  evolution  of  the  energy  carried  by  the  scale 
j  =  6. 

At  the  same  times,  the  same  elements  are  plotted  for  a  calculation  run  using  = 
7  and  a  regular  grid.  Figure  9  reveals  that  the  lack  of  resolution  in  the  strong  gradient 
region  (around  x  =  .5)  induces  numerical  oscillations  of  the  solutions  close  to  those  observed 
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when  a  Gibbs  phenomenon  occurs.  These  oscillations  are  however  localized  in  space  mainly 
thanks  to  the  filtering  operation  performed  by  the  spline  projection  of  the  variables  and  their 
derivatives.  Figure  10  shows  again  the  time  evolution  of  an  approximate  solution  obtained 
with  a  computation  started  with  j m  =  6  (that  is  to  say  26  =  64  degrees  of  freedom  and 
jss  =  5)  and  that  has  been  obtained  using  the  adaptative  procedure.  From  a  time  close  to 
ta  —  .5/tt  when  the  local  gradients  begin  to  be  very  large  the  adaptative  procedure  adapts 
progressively  the  space  to  jMa  =  7  and  ( kltk2 )  =  (22,31)  (that  is  to  say  64  +  20  =  84 
degrees  of  freedom)  and  then  (for  .66/n  <  t  <  IA/tt)  to  jMa  =  8  and  ( k[,k'2 )  =  (54,73) 
(that  is  to  say  64  +  20  +  20  =  104  functions).  The  corresponding  dyadic  grid  is  the  one 
of  Figure  4,  the  smallest  scales  j  =  6  and  j  =  7  being  used  only  when  required.  Finally 
on  Figure  11,  the  evolution  of  the  approximate  solution  using  the  same  time  discretization 
and  a  Fourier  pseudospectral  spatial  discretization  is  plotted.  The  resolution  is  27  =  128. 
First,  one  observes,  in  comparison  to  Figure  9  which  is  obtained  with  the  same  resolution 
(i.e  the  same  number  of  degrees  of  freedom),  that  the  oscillations  are  much  more  spread  in 
the  Fourier  case,  due  to  the  non  localization  of  the  basis.  Secondly,  comparison  with  Figure 
10  shows  that  with  a  maximum  of  104  well  distributed  degrees  of  freedom,  a  much  better 
result  can  be  obtain  from  the  adapted  algorithm.  It  is  known  (see  [Basdevant  et  al.  1986]) 
that  to  improve  the  precision  of  a  Fourier  spectral  method  without  increasing  the  number  of 
degrees  of  freedom,  intricate  coordinate  transformations  are  required  that  imply  moreover 
an  a  priori  knowledge  of  the  solution  shape. 

A  comparison  of  the  four  different  approximations  presented  above  is  made  in  Table  2. 
Tnmx  is  the  time  needed  to  reach  the  maximum  slope  at  x  =  0.5,  and  Sm&x  is  the  value  of 
this  slope.  No  optimization  study  of  the  time  marching  scheme  has  been  made  that  would 
have  improved  the  estimates  of  S^x  and  (see  [Basdevant  et  al.  1986]). 

Table  2:  Comparison  of  some  different  methods. 


Algorithm 

(AB-CN  Time  scheme) 

2t rT^x 
Exact:1.6037 

W2- 

Exact:  152. 005 

Degrees  of 
freedom 

Present  Algorithm,  m  =  6 

3m  -  8,  regular  grid 

1.64 

150.3 

No 

oscillations 

256 

Present  Algorithm,  m  =  6 

3m  =  7,  regular  grid 

1.63 

135.0 

Localized 

oscillations 

128 

Present  Algorithm,  m  =  6 

3u  =  6,;Af„  =  8  adapted  grid 

1.64 

150.3 

No 

oscillations 

<  104 

Fourier  pseudospectral 
n  =  128 

1.62 

134.8 

Spread 

oscillations 

128 
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5.  CONCLUSIONS 


The  numerical  resolution  of  the  ID  periodic  and  regularized  Burgers  equation  has  been 
performed  using  a  new  algorithm.  This  algorithm  is  based  on  the  wavelet  representation  of 
the  space  of  approximation.  In  the  so-called  regular  grid  case,  this  space  is  nothing  more 
than  the  classical  spline  space,  but  in  the  non  regular  grid  case  it  is  more  flexible  than  a 
classical  grid  refinement  generated  space. 

The  special  feature  of  this  algorithm  that  uses  explicitly  the  wavelet  basis  stands  on 
the  fact  that  it  extensively  takes  into  account  the  localization  properties  of  the  wavelets  in 
the  physical  and  Fourier  spaces.  Indeed,  the  exponential  decay  of  the  spline  wavelets  in 
both  spaces  allows,  on  one  hand  handling  the  differential  operators  as  if  they  were  nearly 
diagonal  and,  on  the  other  hand,  working  in  the  physical  space  with  very  flexible  orthogonal 
and  numerically  well  localized  functions.  The  scale  decomposition  provided  by  the  wavelet 
coefficients  and  the  knowledge  of  the  Burgers  operators  axe  used  to  manage  at  each  time 
step  the  space  of  approximation.  The  whole  algorithm  can  then  adapt  itself  to  the  solution. 
Moreover,  in  the  case  of  constant  regular  operators,  advantage  can  be  taken  by  pre-inverting 
the  problem  at  each  scale  once  and  for  all  at  the  beginning  of  the  calculation.  Then,  the 
whole  algorithm  can  be  seen  as  a  time  advancing  pyramidal  algorithm. 

In  this  paper,  the  method  has  been  described  and  a  numerical  application  has  been 
provided  that  uses  spline  wavelets  of  order  6.  Some  mathematical  remarks  have  been  given 
but  the  results  on  numerical  analysis  of  the  method  will  be  published  elsewhere. 

To  date,  we  are  not  sure  that  the  choice  of  the  spline  wavelets  is  optimal  for  efficient 
and  adaptable  resolution  of  partial  differential  equations;  furthermore,  the  existence  of  a 
stronger  connection  between  the  wavelets  and  the  resolved  equation  at  each  scale  has  not 
been  investigated.  However,  it  must  be  recalled  that  the  spline  wavelets  achieve  an  optimal 
combined  localization  in  physical  and  Fourier  spaces,  as  far  as  orthonormal  wavelets  are 
concerned. 

The  limitation  to  a  1-D  problem  and  periodic  boundary  conditions  is  not  connected  to 
the  algorithm  or  the  wavelet  approach  even  if  various  specific  problems  will  have  to  be  faced. 
This  choice  has  been  made  for  clarity  and  simplicity.  Extensions  to  more  realistic  problems 
involving  different  boundary  conditions  and  multidimensions  are  currently  underway. 


18 


REFERENCES 


C.  Basdevant,  M.  Deville,  P.  Haldenwang,  J.  M.  Lacroix,  J.  Ouazzani,  and  R.  Peyret:  “Spec¬ 
tral  and  finite  difference  solution  of  the  Burgers  equation,”  Computers  and  Fluids,  Vol. 
14,  No.  1,  pp.  23-41  (1986). 

G.  Battle:  “A  block  spin  construction  of  ondelettes,  Part  I:  Lemarie  functions,”  Commun. 
Math.  Phys.  110,  601-615  (1987). 

C.  Canuto,  M.  Y.  Hussaini,  A.  Quarteroni,  and  T.  A.  Zang:  “Spectral  methods  in  fluid 

dynamics,”  Springer  Series  in  Computational  Physics  (1987). 

D.  Gottlieb  and  S.  Orszag:  “Numerical  analysis  of  spectral  methods:  Theory  and  applica¬ 

tions,”  SIAM-CBMS  Philadelphia  (1977). 

P.  G.  Lemarie:  “Ondelettes  a  localisation  exponentielles,”  Journ.  de  Math.  Pures  et  Appl. 
1989. 

S.  Mallat:  “Multiresolution  representations  and  wavelets,”  Ph.D.,  MS-CIS-88-68,  Grasp  Lab 
153,  University  of  Pennsylvania,  Philadelphia,  PA  19104  (1988). 

S.  Mallat:  “A  characterization  of  signals  from  the  wavelet  transform  maxima,”  NSF/CBMS 
Conference  on  Wavelets,  Lowell  University  (1990). 

Y.  Meyer:  “Ondelettes  et  Operateurs  1:  Ondelettes,”  Herman  (1990). 

V.  Perrier  and  C.  Basdevant:  “La  decomposition  en  ondelettes  periodiques,  un  outil  pour 
l’analyse  de  champs  inhomogenes.  Theorie  et  algorithmes,”  la  Recherche  Aerospatiale, 
No.  3,  pp.  53-67  (1989). 


19 


Amplitude 


h(x)  Physical  Space 


*) 


OS 


0/1 


0  1 


0.7 


II  I 


mi  •••••••••••' 


o  1  Li _ i _ L 

Pi  i  IS  -  1 


Figure  2. 


•  • 

r>  n 

x 


_i _ i _ i — 

in  is  70  h(ui)  Fourier  Space 


a)  h  function,  for  the  m  =  6  spline  multiresolution  analysis,  in 
physical  and  Fourier  spaces. 

b)  g  function,  for  the  m  =  6  spline  multiresolution  analysis,  in 
physical  and  Fourier  spaces. 
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Time  evolution  of  the  approximated  solution:  m  = 
regular  grid  algorithm.  256  degrees  of  freedom. 


2 


of  the  Approximated  Solution  256  degrees  of  freedom 


Figure  7.  (Continued) 
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of  the  Approximated  Solution  256  degrees  of  freedom 


Figure  9. 


Time  evolution  of  the  approximated  solution  m  —  6,  ju  —  7  regular 
grid  algorithm.  128  degrees  of  freedom. 
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Figure  10.  Time  evolution  of  the  approximated  solution  m  =  6 ,ju  =  6,^.1  = 
7,  adapted  algorithm  .  The  maximum  number  of  degrees  of 
freedom  during  the  whole  computation  is  104. 


33 


NASA 


Report  Documentation  Page 


1  Report  No 

NASA  CR- 187480 
ICASE  Report  No.  90-83 


2  Government  Accession  No 


3.  Recipient's  Catalog  No 


4.  Title  and  Subtitle 

RESOLUTION  OF  THE  ID  REGULARIZED  BURGERS  EQUATION 
USING  A  SPATIAL  WAVELET  APPROXIMATION 


6.  Report  Date 
December  1990 


6  Performing  Organization  Code 


7.  Author(s) 

J.  Liandrat 
Ph.  Tchamitchian 


9  Performing  Organization  Name  and  Address 

Institute  for  Computer  Applications  in  Science 
and  Engineering 

Mail  Stop  132C,  NASA  Langley  Research  Center 
Hampton,  VA  23663-5225 


12.  Sponsoring  Agency  Name  and  Address 

National  Aeronautics  and  Space  Administration 
Langley  Research  Center 
Hampton,  VA  23665-5225 


8.  Performing  Organization  Report  No 

90-83 


10.  Work  Unit  No. 

505-90-21-01 


11.  Contract  or  Grant  No 

NAS1-18605 


13.  Type  of  Report  and  Period  Covered 
Contractor  Report 


14.  Sponsoring  Agency  Code 


15.  Supplementary  Notes 

Langley  Technical  Monitor: 
Richard  W.  Barnwell 


Submitted  to  Annales  Institut 
Henri  Po incard 


Final  Report _ 

16.  Abstract 

The  Burgers  equation  with  a  small  viscosity  term,  initial  and  periodic 
boundary  conditions  is  resolved  using  a  spatial  approximation  constructed  from  an 
orthonormal  basis  of  wavelets. 

The  algorithm  is  directly  derived  from  the  notions  of  multiresolution  analy¬ 
sis  and  tree  algorithms.  Before  the  numerical  algorithm  is  described  these  notions 
are  first  recalled.  The  method  uses  extensively  the  localization  properties  of  the 
wavelets  in  the  physical  and  Fourier  spaces.  Moreover,  we  take  advantage  of  the  fact 
that  the  involved  linear  operators  have  constant  coefficients.  Finally,  the  algo¬ 
rithm  can  be  considered  as  a  time  marching  version  of  the  tree  algorithm. 

The  most  important  point  is  that  an  adaptive  version  of  the  algorithm  exists: 
it  allows  one  to  reduce  in  a  significant  way  the  number  of  degrees  of  freedom  re¬ 
quired  for  a  good  computation  of  the  solution. 

Numerical  results  and  description  of  the  different  elements  of  the  algcrithm 
are  provided  in  combination  with  different  mathematical  comments  on  the  method  and 
some  comparison  with  more  classical  numerical  algorithms. 

17.  Key  Words  (Suggested  by  Author(s))  18.  Distribution  Statement 

wavelet  approximation,  numerical  method,  64  -  Numerical  Analysis 

P.D.E. 


Unclassified  -  Unlimited 


19.  Security  Classil.  lot  this  reportl 

20  Security  Classif.  (of  this  page) 

21  No.  of  pages 

22  Price 

Unclassified 

Unclassified 

36 

A03 

NASA  FORM  1C2S  OCT  88 


NASA-Langley,  1991 


